#As a placebo test this files estimates the effect on past electoral success

rm(list=ls(all=TRUE))
graphics.off()
library(rdd)
library(xtable)
library(dplyr)

setwd() # set working directory
load("data/data_alt_boot.rdata")

#if cases are extremely close they might have wrong signed simscore. Recode these to arbitrarily small simscore
data_alt$margsim[data_alt$margsim <= 0 & data_alt$elected == 1] <-  0.001
data_alt$margsim[data_alt$margsim >= 0 & data_alt$elected == 0] <- -0.001

#subset to 2005 and 2009 data
inc_2005 <- data.frame(subset(data_alt, year == 2005, select = c("muncpr", "candname", "candpart", "elected")))
dat_2009 <- data.frame(subset(data_alt, year == 2009))

dat_2009 <- subset(dat_2009, 
                   marg.los  != 0 & marg.win     != 1 & 
                   openlist  == 1 & abs(margsim) <  1 &
                   splitlist == 0 & 
                   simscore  > 0  & simscore     <  1,
                   select = c("muncpr", "candname", "candpart", "margsim"))

#merge data frames
data_plac <- left_join(dat_2009, inc_2005, by = c("muncpr", "candname", "candpart"))

#create variable for running in t-1
data_plac$ran     <- 1 - is.na(data_plac$elected)
#recode variable for being elected in t-1
data_plac$elected[is.na(data_plac$elected)] <- 0


##Estimate and plot RD 
##Incumbent 
rdest    <- RDestimate(elected ~ margsim, data = data_plac)
rdest.cl <- RDestimate(elected ~ margsim, data = data_plac, cluster = data_plac$muncpr)

graphics.off()
quartz(w=8, h =4)
par(mfrow=c(1,2), 
    mar=c(4,4,2,1),
    cex = 1.1)

plot(rdest, range = c(-.8,.8))
title( ylab = "Pr(elected in 2005)", xlab = "Votes relative to threshold")
#add rugplot to illustrate point density over running var
rugz1 <- quantile(data_plac$margsim, probs = 0:100/100, na.rm = TRUE)
rug(rugz1)
abline(v=0, lty = 2, col = "grey50")

#plot for Ran
rdest1    <- RDestimate(ran ~ margsim, data = data_plac)
rdest1.cl <- RDestimate(ran ~ margsim, data = data_plac, cluster = data_plac$muncpr)
plot(rdest1, range = c(-.8,.8))
title( ylab = "Pr(ran in 2005)", xlab = "Votes relative to threshold")
rug(rugz1)
abline(v=0, lty = 2, col = "grey50")

##Summarize jumps in table

out1 <- cbind(round(rdest$est,2),
              round(rdest$se,2),
              round(rdest.cl$se,2),
              round(rdest$bw,2),
              rdest$obs)
colnames(out1) <- c("Estimate", "Std. Error", "Rob. Std. Error", "Bandwidth", "Observations")

out2 <- cbind(round(rdest1$est,2),
              round(rdest1$se,2),
              round(rdest1.cl$se,2),
              round(rdest1$bw,2),
              rdest1$obs)
colnames(out2) <- c("Estimate", "Std. Error", "Rob. Std. Error", "Bandwidth", "Observations")

xtable(cbind(t(out1),t(out2)))
